% solves set of 9 differential equations used during cardiac cycle

function dy = eject(t,y,flag,R,C,L,Z0,time,E_t_LV,dEdt_LV,Rao,Rmv,VdeadV,Pvenous);

dy=zeros(3,1);

% y1 = VLV
% y2 = Pr2
% y3 = Pao

% Flows and volumes are calculated directly from the pressures and resistances

% p1 is value of LV E(t) at timestep
% p3 is value of LV dE(t) at timestep

p1=interp1([time' max(time+0.001)]',[E_t_LV' E_t_LV(length(E_t_LV))]',t);
p3=interp1([time' max(time+0.001)]',[dEdt_LV' dEdt_LV(length(dEdt_LV))]',t);

dummyRmv=Rmv(sign(p1*(y(1)-VdeadV)-Pvenous)+2);
dummyRao=Rao(sign(y(3)-p1*(y(1)-VdeadV))+2);

% Qmv=(Pvenous-PLV)/dummyRmv;
% Qao=(PLV-y(4))/dummyRao;
% Qpv=(Pvenous-PLA)/Rpv;

% PLV=p1*(y(1)-VdeadV);
% PLA=p2*(y(2)-VdeadA);

dy(1)=(Pvenous-p1*(y(1)-VdeadV))/dummyRmv-(p1*(y(1)-VdeadV)-y(3))/dummyRao;;
dy(2)=((p1*(y(1)-VdeadV)-y(3))/dummyRao)/C-y(2)/(R*C);																		%dPr2/dt = Qao/C - Pr2/(RC)
dy(3)=(Z0*(p1*dy(1)+p3*(y(1)-VdeadV))/dummyRao+dy(2)-Z0*y(3)/L+Z0*y(2)/L)/(1+Z0/dummyRao); 		%dPao/dt = Z0 dQao/dt + dPr2/dt - Z0*Pao/L + Z0*Pr2/L

% initial procedure: perfectly closing mitral and aortic valve
% forward resistance = finite; backward resistance is infinite
% data are stored in resistance-matrix
% first element = forward resistance; 3rd element is backward resistance
% values are switched based on sign of PLV-PLA and PAO-PLV 
% if sign = -1: forward flow
% if sign = 1: backward flow
% resistance = resistance(sign(**)+2)
% for forward flow, the 1st element is used; for backward flow the 3rd element is used
